####Overview####
# File name:    replication_jepp
# Last Change:  27/02/2024
# Name:         David Hope/Julian Limberg/Nina Weber
# Purpose:      Replication Data for "The ICT Revolution and Preferences for Taxing Top Earners" (JEPP)

####1. System setups#####
#Load Packages
# install packages from CRAN
p_needed <- c("readstata13", "tidyverse", "texreg","varhandle","ggthemes","grid","cowplot",
              "foreign","magrittr","haven","xtable","ggrepel","ggplot2","ggpubr", "fixest",
              "Matrix","gridExtra","grid","broom", "mediation",
              "plotly","readxl","ggdag","dagitty","xlsx","scales","estimatr","miceadds")
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}
lapply(p_needed, require, character.only = TRUE)

#Set Seed
set.seed(2789)

#Set working directory to current folder
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


load("jepp_data.RData")



spec_red <- spec %>%
  filter(treatment != 0)

#Drop Worker who claims to be younger than 18

workers_red <- workers %>%
  filter(treatment != 0,
         d1 >= 18)

workers2_win <- workers%>%
  mutate(m2_change = (m2b-m2a)/m2a) %>%
  filter(position == 1)

workers2_lose <- workers%>%
  mutate(m2_change = (m2b-m2a)/m2a) %>%
  filter(position == 0)

workers2_win_red <- workers2_win %>%
  filter(treatment != 0)

workers2_lose_red <- workers2_lose %>%
  filter(treatment != 0)



####2. Spectators -- Distribution Points####

spec_luck <- spec %>%
  filter(treatment == 0) %>%
  pivot_longer(c(b7_1,b7_2,b7_3,b7_4), names_to = "Traits", values_to = "Characteristics") %>%
  dplyr::mutate(Traits = case_when(
    Traits == "b7_1" ~ "Inherited Intelligence",
    Traits == "b7_2" ~ "Education",
    Traits == "b7_3" ~ "Luck",
    Traits == "b7_4" ~ "Effort"
  ),
  Treatment = "Luck") %>%
  dplyr::group_by(Traits) %>%
  summarise(mean_Characteristics = mean(Characteristics,na.rm = T)) %>%
  mutate(Treatment = "Luck")


spec_luck$Traits <- factor(spec_luck$Traits, levels = c("Luck", "Effort", "Education", "Inherited Intelligence"))

spec_routine <- spec %>%
  filter(treatment == 1) %>%
  pivot_longer(c(b7_1,b7_2,b7_3,b7_4), names_to = "Traits", values_to = "Characteristics") %>%
  dplyr::mutate(Traits = case_when(
    Traits == "b7_1" ~ "Inherited Intelligence",
    Traits == "b7_2" ~ "Education",
    Traits == "b7_3" ~ "Luck",
    Traits == "b7_4" ~ "Effort"
  ),
  Treatment = "Routine Work") %>%
  dplyr::group_by(Traits) %>%
  summarise(mean_Characteristics = mean(Characteristics,na.rm = T)) %>%
  mutate(Treatment = "Routine Work")

spec_routine$Traits <- factor(spec_routine$Traits, levels = c("Luck", "Effort", "Education", "Inherited Intelligence"))


spec_complex <- spec %>%
  filter(treatment == 2) %>%
  pivot_longer(c(b7_1,b7_2,b7_3,b7_4), names_to = "Traits", values_to = "Characteristics") %>%
  dplyr::mutate(Traits = case_when(
    Traits == "b7_1" ~ "Inherited Intelligence",
    Traits == "b7_2" ~ "Education",
    Traits == "b7_3" ~ "Luck",
    Traits == "b7_4" ~ "Effort"
  ),
  Treatment = "Complex Work") %>%
  dplyr::group_by(Traits) %>%
  summarise(mean_Characteristics = mean(Characteristics,na.rm = T)) %>%
  mutate(Treatment = "Complex Work")

spec_complex$Traits <- factor(spec_complex$Traits, levels = c("Luck", "Effort", "Education", "Inherited Intelligence"))



luck_fig <- ggplot(spec_luck, aes(x= Traits, y= mean_Characteristics)) +
  geom_bar(stat="identity") +
  ylab("Point Allocation") +
  ggtitle("Luck") +
  ylim(0,100) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank(),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(size=20))


spec_routine <- spec %>%
  filter(treatment == 1) %>%
  pivot_longer(c(b7_1,b7_2,b7_3,b7_4), names_to = "Traits", values_to = "Characteristics") %>%
  dplyr::mutate(Traits = case_when(
    Traits == "b7_1" ~ "Inherited Intelligence",
    Traits == "b7_2" ~ "Education",
    Traits == "b7_3" ~ "Luck",
    Traits == "b7_4" ~ "Effort"
  )) %>%
  dplyr::group_by(Traits) %>%
  summarise(mean_Characteristics = mean(Characteristics,na.rm = T))

spec_routine$Traits <- factor(spec_routine$Traits, levels = c("Luck", "Effort", "Education", "Inherited Intelligence"))


routine_fig <- ggplot(spec_routine, aes(x= Traits, y= mean_Characteristics)) +
  geom_bar(stat="identity") +
  ggtitle("Routine Work") +
  ylab("Point Allocation") +
  ylim(0,100) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank(),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(size=20))
spec_complex <- spec %>%
  filter(treatment == 2) %>%
  pivot_longer(c(b7_1,b7_2,b7_3,b7_4), names_to = "Traits", values_to = "Characteristics") %>%
  dplyr::mutate(Traits = case_when(
    Traits == "b7_1" ~ "Inherited Intelligence",
    Traits == "b7_2" ~ "Education",
    Traits == "b7_3" ~ "Luck",
    Traits == "b7_4" ~ "Effort"
  )) %>%
  dplyr::group_by(Traits) %>%
  summarise(mean_Characteristics = mean(Characteristics,na.rm = T))

spec_complex$Traits <- factor(spec_complex$Traits, levels = c("Luck", "Effort", "Education", "Inherited Intelligence"))


complex_fig <- ggplot(spec_complex, aes(x= Traits, y= mean_Characteristics)) +
  geom_bar(stat="identity") +
  ggtitle("Complex Work") +
  ylab("Point Allocation") +
  ylim(0,100) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank(),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        text = element_text(size=20))


combine_descriptive <- ggarrange(luck_fig, routine_fig, complex_fig, nrow = 1)

####3. Main Analysis -- Both####




spec_min <- spec %>%
  dplyr::select(rd, treatment) %>%
  mutate(Group = "Spectators",
         treatment = case_when(
           treatment == 0 ~ "Luck",
           treatment == 1 ~ "Routine Work",
           treatment == 2 ~ "Complex Work"
         ))

workers_min <- workers %>%
  dplyr::select(rd, treatment) %>%
  mutate(Group = "Workers",
         treatment = case_when(
           treatment == 0 ~ "Luck",
           treatment == 1 ~ "Routine Work",
           treatment == 2 ~ "Complex Work"
         ))


both_summary <- rbind(spec_min, workers_min)

both_summary$treatment <- factor(
  both_summary$treatment, 
  levels = c("Luck", "Routine Work","Complex Work"))




mean <- both_summary %>% 
  group_by(Group, treatment)%>%
  dplyr::summarise(mean_rd = mean(rd, na.rm= T))

sd <- both_summary %>% 
  group_by(Group, treatment)%>%
  dplyr::summarise(sd_rd = sd(rd, na.rm= T))

cbind(mean,sd) %>%
  mutate((sd_rd/mean_rd)*100)

predicted_rate <- mean %>%
  ggplot(., aes(x = factor(treatment), y=mean_rd, fill = Group)) +
  geom_bar(position = "dodge",stat="identity") + 
  theme_bw()+ 
  labs(y = "Mean Top Tax Rate",
       x = "Treatment Group") + 
  scale_fill_grey() +
  theme(text = element_text(size=20),
        legend.text=element_text(size=15),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.title=element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        legend.position = c(.8,.9))

#Analysis


summary(model_1 <- lm_robust(rd ~ factor(treatment), spec, clusters = prolific))
summary(model_2 <- lm_robust(rd ~ factor(treatment), spec_red, clusters = prolific))


summary(model_3 <- lm_robust(data=workers, rd ~ factor(treatment)))
summary(model_4 <- lm_robust(data=workers_red, rd ~ factor(treatment)))


est_all_1 <- tidy(model_1, conf.int = TRUE) %>%
  mutate(model = "Reference Group: Luck",
         Group = "Spectators") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))
est_all_2 <- tidy(model_2, conf.int = TRUE) %>%
  mutate(model = "Reference Group: Routine Work",
         Group = "Spectators") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))

est_all_3 <- tidy(model_3, conf.int = TRUE) %>%
  mutate(model = "Reference Group: Luck",
         Group = "Workers") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))
est_all_4 <- tidy(model_4, conf.int = TRUE) %>%
  mutate(model = "Reference Group: Routine Work",
         Group = "Workers") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))


screenreg(l = list(model_1, model_2, model_3, model_4), include.ci = FALSE,
       include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
       label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
       custom.coef.map = list("factor(treatment)1" = "Routine Work",
                              "factor(treatment)2"= "Complex Work"))



all_treat <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
  dplyr::mutate(term = case_when(
    term == "factor(treatment)1" ~ "Routine Work",
    term == "factor(treatment)2" ~ "Complex Work"
  )) %>%
  mutate(signif_stars = case_when(
    p.value < 0.05 & p.value > 0.01~ "*",
    p.value < 0.01 & p.value > 0.001 ~ "**",
    p.value < 0.001 ~ "***",
    T ~ ""))

est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
  mutate(model = "Reference Group: Luck",
         Group = "Spectators") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))
est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
  mutate(model = "Reference Group: Routine Work",
         Group = "Spectators") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))

est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
  mutate(model = "Reference Group: Luck",
         Group = "Workers") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))
est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
  mutate(model = "Reference Group: Routine Work",
         Group = "Workers") %>%
  filter(term %in% c("factor(treatment)1","factor(treatment)2"))


all_treat_2 <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
  dplyr::mutate(term = case_when(
    term == "factor(treatment)1" ~ "Routine Work",
    term == "factor(treatment)2" ~ "Complex Work"
  ))


all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model", "Group"))



all_treat$term <- factor(
  all_treat$term, 
  levels = c("Complex Work","Routine Work"))



all_plot <- all_treat %>%
  filter(model == "Reference Group: Luck") %>%
  ggplot(., aes(x=term, y=estimate, color = Group))  +
  geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+ 
  geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+ 
  geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
  theme_bw()+ 
  geom_point(size=3, position = position_dodge(width = 0.8)) +
  labs(y = "Treatment Effect on Tax Rate",
       x = "") + 
  scale_colour_grey() +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(breaks=seq(-40,10,10)) + 
  ylim(-40,10) +
  coord_flip() +
  theme(text = element_text(size=20),
        axis.title.x=element_blank(),
        legend.text=element_text(size=15),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.position = "none") +
  facet_wrap(.~model)





reduced_plot <- all_treat %>%
  filter(model == "Reference Group: Routine Work") %>%
  ggplot(., aes(x=term, y=estimate, color = Group))  +
  geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
  geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
  geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
  theme_bw()+ 
  geom_point(size=3, position = position_dodge(width = 0.8)) +
  labs(y = "Treatment Effect on Tax Rate",
       x = "") +
  scale_colour_grey() +
  geom_hline(yintercept=0, linetype = "dashed") +
  scale_y_continuous(breaks=seq(-40,10,10)) + 
  ylim(-40,10) +
  coord_flip() +
  theme(text = element_text(size=20),
        legend.text=element_text(size=15),
        plot.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank(),
        legend.title=element_blank(),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"),
        legend.position = "bottom")  +
  facet_wrap(.~model)


combine_descriptive <- ggarrange(all_plot, reduced_plot, nrow = 2, heights=c(3,4))


####4. Interaction####

#Male

spec_red <- spec_red %>%
  mutate(gender = case_when(
    d2 == "Male" ~ 1,
    d2 == "Female" ~ 0,
    T ~ NA_real_
  ))

workers_red <- workers_red %>%
  mutate(gender = case_when(
    d2 == "Male" ~ 1,
    d2 == "Female" ~ 0,
    T ~ NA_real_
  ))





summary(model_1 <- lm_robust(rd ~ factor(treatment)* gender, spec_red, clusters = prolific))


summary(model_2 <- lm_robust(data=workers_red, rd ~ factor(treatment)* gender))


#Age
summary(model_3 <- lm_robust(rd ~ factor(treatment)* d1, spec_red, clusters = prolific))


summary(model_4 <- lm_robust(data=workers_red, rd ~ factor(treatment)* d1))

#Republican

spec_red <- spec_red %>%
  mutate(party = case_when(
    d9 == "Republican Party" ~ 1,
    d9 == "Democratic Party" ~ 0,
    T ~ NA_real_
  ))
  
  workers_red <- workers_red %>%
    mutate(party = case_when(
      d9 == "Republican Party" ~ 1,
      d9 == "Democratic Party" ~ 0,
      T ~ NA_real_
    ))
  
  
  summary(model_5 <- lm_robust(rd ~ factor(treatment)* party, spec_red, clusters = prolific))
  
  
  summary(model_6 <- lm_robust(data=workers_red, rd ~ factor(treatment)* party))
  
  #Income%>%
  
  
  spec_red <- spec_red %>%
    mutate(income = case_when(
      d5 == "Under $10,000" ~ 5000,
      d5 == "$10,000 - $20,000" ~ 15000,
      d5 == "$20,001 - $30,000" ~ 25000,
      d5 == "$30,001 - $40,000" ~ 35000,
      d5 == "$40,001 - $50,000" ~ 45000,
      d5 == "$50,001 - $60,000" ~ 55000,
      d5 == "$60,001 - $80,000" ~ 70000,
      d5 == "$80,001 - $100,000" ~ 90000,
      d5 == "$100,001 - $150,000" ~ 125000,
      d5 == "$150,001 - $200,000" ~ 175000,
      d5 == "$200,001 - $350,000" ~ 275000,
      d5 == "$350,001 - $500,000" ~ 425000,
      T ~ NA_real_))
  
  workers_red <- workers_red %>%
    mutate(income = case_when(
      d5 == "Under $10,000" ~ 5000,
      d5 == "$10,000 - $20,000" ~ 15000,
      d5 == "$20,001 - $30,000" ~ 25000,
      d5 == "$30,001 - $40,000" ~ 35000,
      d5 == "$40,001 - $50,000" ~ 45000,
      d5 == "$50,001 - $60,000" ~ 55000,
      d5 == "$60,001 - $80,000" ~ 70000,
      d5 == "$80,001 - $100,000" ~ 90000,
      d5 == "$100,001 - $150,000" ~ 125000,
      d5 == "$150,001 - $200,000" ~ 175000,
      d5 == "$200,001 - $350,000" ~ 275000,
      d5 == "$350,001 - $500,000" ~ 425000,
      T ~ NA_real_))
  
  
  summary(model_7 <- lm_robust(rd ~ factor(treatment)* income, spec_red, clusters = prolific))
  
  
  summary(model_8 <- lm_robust(data=workers_red, rd ~ factor(treatment)* income))
  
  
  #Education
  
  
  spec_red <- spec_red %>%
    mutate(college = case_when(
      d4 == "2-year college degree" ~1,
      d4 == "4-year college degree" ~1,
      d4 == "Doctoral degree" ~1,
      d4 == "Professional degree (JD, MD, MBA)" ~1,
      is.na(d4) ~ NA_real_,
      T ~0))
  
  summary(model_9 <- lm_robust(rd ~ factor(treatment)* college, spec_red, clusters = prolific))
  
  
  texreg(l = list(model_1, model_3, model_5, model_7, model_9), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_beliefs", caption = "Regression Results for Interaction Effects", caption.above = TRUE, 
         custom.coef.map = list("factor(treatment)2"= "Complex Work",
                                "gender" = "Male",
                                "factor(treatment)2:gender" = "Complex Work * Male",
                                "d1" = "Age",
                                "factor(treatment)2:d1" = "Complex Work * Age",
                                "party" = "Republican",
                                "factor(treatment)2:party" = "Complex Work * Republican",
                                "income" = "Income",
                                "factor(treatment)2:income" = "Complex Work * Income",
                                "college" = "College Degree",
                                "factor(treatment)2:college" = "Complex Work * College Degree"))
  
  
  
####5. Spectators -- Beliefs####
  summary(model_1 <- lm_robust(data=spec, b6 ~ factor(treatment), clusters = prolific))
  summary(model_2 <- lm_robust(data=spec_red, b6 ~ factor(treatment), clusters = prolific))
  summary(model_3 <- lm_robust(data=spec, b4 ~ factor(treatment), clusters = prolific))
  summary(model_4 <- lm_robust(data=spec_red, b4 ~ factor(treatment), clusters = prolific))
  summary(model_5 <- lm_robust(data=spec, b5 ~ factor(treatment), clusters = prolific))
  summary(model_6 <- lm_robust(data=spec_red, b5 ~ factor(treatment), clusters = prolific))
  summary(model_7 <- lm_robust(data=spec, b2 ~ factor(treatment), clusters = prolific))
  summary(model_8 <- lm_robust(data=spec_red, b2 ~ factor(treatment), clusters = prolific))
  summary(model_9 <- lm_robust(data=spec, b3 ~ factor(treatment), clusters = prolific))
  summary(model_10 <- lm_robust(data=spec_red, b3 ~ factor(treatment), clusters = prolific))

  #Only Work (Routine vs Complex)
  
  est_all_1 <- tidy(model_2, conf.int = TRUE) %>%
    mutate(model = "Luck Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_4, conf.int = TRUE) %>%
    mutate(model = "Effort Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_3 <- tidy(model_6, conf.int = TRUE) %>%
    mutate(model = "Skill Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_8, conf.int = TRUE) %>%
    mutate(model = "Fairness Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_5 <- tidy(model_10, conf.int = TRUE) %>%
    mutate(model = "Deservingness Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  all_treat <- rbind(est_all_1, est_all_2,est_all_3, est_all_4, est_all_5) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ),
    signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  est_all_6 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Luck Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_7 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Effort Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_8 <- tidy(model_6, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Skill Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_9 <- tidy(model_8, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Fairness Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_10 <- tidy(model_10, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Deservingness Belief") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  all_treat_2 <- rbind(est_all_6, est_all_7,est_all_8, est_all_9, est_all_10) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  all_treat_2$term <- factor(
    all_treat_2$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model"))
  
  
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  plot_1 <- all_treat %>%
    filter(model == "Luck Belief") %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Treatment Effect on Luck Belief",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank()) +
    facet_wrap(.~model)
  
  plot_2 <- all_treat %>%
    filter(model == "Effort Belief") %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Treatment Effect on Effort Belief",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank()) +
    facet_wrap(.~model)
  
  plot_3 <- all_treat %>%
    filter(model == "Skill Belief") %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Treatment Effect on Skill Belief",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank()) +
    facet_wrap(.~model)
  
  plot_4 <- all_treat %>%
    filter(model == "Fairness Belief") %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Treatment Effect on Fairness Belief",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank()) +
    facet_wrap(.~model)
  
  plot_5 <- all_treat %>%
    filter(model == "Deservingness Belief") %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Treatment Effect on Deservingness Belief",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank()) +
    facet_wrap(.~model)
  
  combine_descriptive <- ggarrange(plot_1, plot_2, plot_3, plot_4, plot_5,nrow = 5)
  

####6. Causal Mediation####
  
  
  
  spec_red_all <- spec_red %>%
    filter(
      !(is.na(b2)),  
      !(is.na(b3)),
      !(is.na(b4)),
      !(is.na(b5)),
      !(is.na(rd))
    ) %>%
    mutate(factor_treatment = factor(treatment)) %>%
    mutate(prolific_cluster = factor(prolific)) %>%
    data.frame()
  
  
  ###Deservigness
  #Treat -> Effort -> Deserving
  fitM1 <- lm(b4 ~ factor_treatment , data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY2 <- lm(b3 ~ factor_treatment + b4 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  cluster_prol <- spec_red_all$prolific %>% c()
  
  fitMed1 <- mediate(fitM1, fitY2, treat="factor_treatment", mediator="b4", cluster = cluster_prol)
  
  summary(fitMed1)
  
  effort_1 <- tidy(fitMed1, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(7.567),
                               p.value = c(0),
                               conf.low = c(5.030),
                               conf.high = c(10)
  )
  
  
  effort_des <- rbind(effort_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_1 <- effort_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Perception of Deservingness (Effort as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  #Treat -> Skill -> Deserving
  fitM3 <- lm(b5 ~ factor_treatment, data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY4 <- lm(b3 ~ factor_treatment + b5 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  fitMed2 <- mediate(fitM3, fitY4, treat="factor_treatment", mediator="b5", cluster = cluster_prol)
  
  summary(fitMed2)
  
  skill_1 <- tidy(fitMed2, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(7.649),
                               p.value = c(2e-16),
                               conf.low = c(5.259),
                               conf.high = c(10.23)
  )
  
  
  skill_des <- rbind(skill_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_2 <- skill_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Perception of Deservingness (Skill as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  
  #Treat -> Deserving -> Red
  fitM5 <- lm(b3 ~ factor_treatment , data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY6 <- lm(rd ~ factor_treatment + b3 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  fitMed3 <- mediate(fitM5, fitY6, treat="factor_treatment", mediator="b3", cluster = cluster_prol)
  
  summary(fitMed3)
  
  red_1 <- tidy(fitMed3, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(-5.290),
                               p.value = c(2e-16),
                               conf.low = c(-8.121),
                               conf.high = c(-2.52)
  )
  
  
  red_des <- rbind(red_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_3 <- red_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Tax Preference (Deservingness as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  combine_med <- ggarrange(fig_1, fig_2, fig_3, nrow = 3, heights=c(4,4))

  ###Fairness
  #Treat -> Effort -> Fair
  fitM1 <- lm(b4 ~ factor_treatment , data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY2 <- lm(b2 ~ factor_treatment + b4 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  cluster_prol <- spec_red_all$prolific %>% c()
  
  fitMed1 <- mediate(fitM1, fitY2, treat="factor_treatment", mediator="b4", cluster = cluster_prol)
  
  summary(fitMed1)
  
  effort_1 <- tidy(fitMed1, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(3.8719),
                               p.value = c(0),
                               conf.low = c(1.1099),
                               conf.high = c(6.58)
  )
  
  
  effort_des <- rbind(effort_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_1 <- effort_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Perception of Fairness (Effort as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  #Treat -> Skill -> Fair
  fitM3 <- lm(b5 ~ factor_treatment, data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY4 <- lm(b2 ~ factor_treatment + b5 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  fitMed2 <- mediate(fitM3, fitY4, treat="factor_treatment", mediator="b5", cluster = cluster_prol)
  
  summary(fitMed2)
  
  skill_1 <- tidy(fitMed2, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(3.80),
                               p.value = c(2e-16),
                               conf.low = c(1.04),
                               conf.high = c(6.63)
  )
  
  
  skill_des <- rbind(skill_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_2 <- skill_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Perception of Fairness (Skill as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  
  #Treat -> Fair -> Red
  fitM5 <- lm(b2 ~ factor_treatment , data=spec_red_all) #IV on M; Hours since dawn predicting coffee consumption
  fitY6 <- lm(rd ~ factor_treatment + b2 , data=spec_red_all) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
  
  fitMed3 <- mediate(fitM5, fitY6, treat="factor_treatment", mediator="b2", cluster = cluster_prol)
  
  summary(fitMed3)
  
  red_1 <- tidy(fitMed3, conf.int = TRUE) %>%
    dplyr::filter(term == "acme_1" | term == "ade_1")%>% dplyr::select(!std.error)
  
  df_full_effect <- data.frame(term = c("Total Effect"),
                               estimate = c(-5.383),
                               p.value = c(2e-16),
                               conf.low = c(-8.198),
                               conf.high = c(-2.40)
  )
  
  
  red_des <- rbind(red_1, df_full_effect) %>%
    mutate(term = case_when(
      term == "acme_1" ~ "ACME",
      term == "ade_1" ~ "ADE",
      term == "Total Effect" ~ "Total Effect",
    ),signif = case_when(
      p.value < 0.05 ~ 1,
      T ~ 0),
    signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  fig_3 <- red_des %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low, ymax=conf.high), width=0, size = 1.2)+
    theme_bw()+ 
    geom_text(position = position_dodge(width = 0.5), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Effect Treatment on Tax Preference (Fairness as Mediator)",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-11,11) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  combine_med <- ggarrange(fig_1, fig_2, fig_3, nrow = 3, heights=c(4,4))
  
  
####7. Specific Mechanisms####
  
  summary(model_1 <- lm_robust(data=spec_red, m1 ~ factor(treatment), clusters = prolific))
  
  
  est_all_1 <- tidy(model_1, conf.int = TRUE) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  est_all_2 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(est_all_1, est_all_2, by = c("term", "estimate", "std.error", "statistic", "p.value"))
  
  fig_1 <- all_treat %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Minimum Bonus Treatment",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  summary(model_2 <- lm_robust(data=spec_red, m2_change ~ factor(treatment), clusters = prolific))
  
  est_all_1 <- tidy(model_2, conf.int = TRUE) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(est_all_1, est_all_2, by = c("term", "estimate", "std.error", "statistic", "p.value"))
  
  fig_2 <- all_treat %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Percentage Difference Between Prize Scenarios",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  
  summary(model_3 <- lm_robust(data=spec_red, m3 ~ factor(treatment), clusters = prolific))
  
  est_all_1 <- tidy(model_3, conf.int = TRUE) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  est_all_2 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2")) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(est_all_1, est_all_2, by = c("term", "estimate", "std.error", "statistic", "p.value"))
  
  fig_3 <- all_treat %>%
    ggplot(., aes(x=term, y=estimate))  +
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.5), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.5)) +
    labs(y = "Workers Above Performance Threshold",
         x = "") + 
    geom_hline(yintercept=0, linetype = "dashed") +
    ylim(-25,25) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank())
  
  combine_descriptive <- ggarrange(fig_1, fig_2, fig_3, nrow = 3, heights=c(1,1,1))
  
 screenreg(l = list(model_1, model_2, model_3), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "mechanisms", caption = "Regression Results for Specific Mechanisms", caption.above = TRUE, 
         custom.coef.map = list("factor(treatment)2"= "Complex Work"))
  
  
####8. Drop quickest 10%####
  quantile(spec$durationinseconds, probs = seq(0, 1, 0.1), na.rm = FALSE)
  quantile(workers$durationinseconds, probs = seq(0, 1, 0.1), na.rm = FALSE)
  
  spec_excl <- spec %>%
    filter(durationinseconds >= 609)
  workers_excl <- workers %>%
    filter(durationinseconds >= 222)
  
  
  spec_red1 <- spec_excl %>%
    filter(treatment != 0)
  workers_red1 <- workers %>%
    filter(treatment != 0)
  
  
  
  summary(model_1 <- lm_robust(rd ~ factor(treatment), spec_excl, clusters = prolific))
  summary(model_2 <- lm_robust(rd ~ factor(treatment), spec_red1, clusters = prolific))
  
  
  summary(model_3 <- lm_robust(data=workers_excl, rd ~ factor(treatment)))
  summary(model_4 <- lm_robust(data=workers_red1, rd ~ factor(treatment)))
  
  
  est_all_1 <- tidy(model_1, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
 screenreg(l = list(model_1, model_2, model_3, model_4), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
         custom.coef.map = list("factor(treatment)1" = "Routine Work",
                                "factor(treatment)2"= "Complex Work"))
  
  
  
  all_treat <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    )) %>%
    mutate(signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  all_treat_2 <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model", "Group"))
  
  
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  
  all_plot <- all_treat %>%
    filter(model == "Reference Group: Luck") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+ 
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+ 
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") + 
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-40,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          axis.title.x=element_blank(),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = "none") +
    facet_wrap(.~model)
  
  
  
  
  
  reduced_plot <- all_treat %>%
    filter(model == "Reference Group: Routine Work") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") +
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-40,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.title=element_blank(),
          legend.background = element_blank(),
          legend.box.background = element_rect(colour = "black"),
          legend.position = "bottom")  +
    facet_wrap(.~model)
  
  
  combine_descriptive <- ggarrange(all_plot, reduced_plot, nrow = 2, heights=c(3,4))

  
####9. Kick Out Those who didn't say 100% luck####
  spec_not_100 <- spec %>%
    filter(treatment == 0) %>%
    filter(b7_3 != 100) %>% 
    select(prolific_pid)
  
  workers_not_100 <- workers %>%
    filter(treatment == 0) %>%
    filter(b7_3 != 100) %>% 
    select(prolific_pid)
  
  spec_100 <- spec %>%
    filter(!(prolific_pid %in% spec_not_100$prolific_pid))
  workers_100 <- workers %>%
    filter(!(prolific_pid %in% workers_not_100$prolific_pid))  
  
  
  spec_red2 <- spec_100 %>%
    filter(treatment != 0)
  workers_red2 <- workers_100 %>%
    filter(treatment != 0)
  
  
  
  summary(model_1 <- lm_robust(rd ~ factor(treatment), spec_100, clusters = prolific))
  summary(model_2 <- lm_robust(rd ~ factor(treatment), spec_red2, clusters = prolific))
  
  
  summary(model_3 <- lm_robust(data=workers_100, rd ~ factor(treatment)))
  summary(model_4 <- lm_robust(data=workers_red2, rd ~ factor(treatment)))
  
  
  est_all_1 <- tidy(model_1, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
screenreg(l = list(model_1, model_2, model_3, model_4), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
         custom.coef.map = list("factor(treatment)1" = "Routine Work",
                                "factor(treatment)2"= "Complex Work"))
  
  
  
  all_treat <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    )) %>%
    mutate(signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  all_treat_2 <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model", "Group"))
  
  
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  
  all_plot <- all_treat %>%
    filter(model == "Reference Group: Luck") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+ 
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+ 
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") + 
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-50,10,10)) + 
    ylim(-50,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          axis.title.x=element_blank(),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = "none") +
    facet_wrap(.~model)
  
  
  
  
  
  reduced_plot <- all_treat %>%
    filter(model == "Reference Group: Routine Work") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") +
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-50,10,10)) + 
    ylim(-50,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.title=element_blank(),
          legend.background = element_blank(),
          legend.box.background = element_rect(colour = "black"),
          legend.position = "bottom")  +
    facet_wrap(.~model)
  
  
  combine_descriptive <- ggarrange(all_plot, reduced_plot, nrow = 2, heights=c(3,4))
  
  
  
  

  
####11. Robust: Losers vs winners####
  
  
  
  summary(model_1 <- lm_robust(rd ~ factor(treatment), workers2_win))
  summary(model_2 <- lm_robust(rd ~ factor(treatment), workers2_win_red))
  
  
  summary(model_3 <- lm_robust(rd ~ factor(treatment), workers2_lose))
  summary(model_4 <- lm_robust(rd ~ factor(treatment), workers2_lose_red))
  
  
  est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.95) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.95) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.95) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Non-Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.95) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Non-Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  screenreg(l = list(model_1, model_2, model_3, model_4), include.ci = FALSE,
            include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
            booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
            label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
            custom.coef.map = list("factor(treatment)1" = "Routine Work",
                                   "factor(treatment)2"= "Complex Work"))
  
  
  
  all_treat <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    )) %>%
    mutate(signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Non-Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Non-Top Earners") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  all_treat_2 <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model", "Group"))
  
  
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  
  all_plot <- all_treat %>%
    filter(model == "Reference Group: Luck") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+ 
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+ 
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") + 
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-30,20) +
    coord_flip() +
    theme(text = element_text(size=20),
          axis.title.x=element_blank(),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = "none") +
    facet_wrap(.~model)
  
  
  
  
  
  reduced_plot <- all_treat %>%
    filter(model == "Reference Group: Routine Work") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") +
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-30,20) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.title=element_blank(),
          legend.background = element_blank(),
          legend.box.background = element_rect(colour = "black"),
          legend.position = "bottom")  +
    facet_wrap(.~model)
  
  
  combine_descriptive <- ggarrange(all_plot, reduced_plot, nrow = 2, heights=c(3,4))
  
  
####12. Robust: Non-Clustered SEs Main####
  
  
  
  summary(model_1 <- lm_robust(rd ~ factor(treatment), spec, se_type = "stata"))
  summary(model_2 <- lm_robust(rd ~ factor(treatment), spec_red, se_type = "stata"))
  
  
  summary(model_3 <- lm_robust(data=workers, rd ~ factor(treatment), se_type = "stata"))
  summary(model_4 <- lm_robust(data=workers_red, rd ~ factor(treatment), se_type = "stata"))
  
  
  est_all_1 <- tidy(model_1, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
screenreg(l = list(model_1, model_2, model_3, model_4), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
         custom.coef.map = list("factor(treatment)1" = "Routine Work",
                                "factor(treatment)2"= "Complex Work"))
  
  
  
  all_treat <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    )) %>%
    mutate(signif_stars = case_when(
      p.value < 0.05 & p.value > 0.01~ "*",
      p.value < 0.01 & p.value > 0.001 ~ "**",
      p.value < 0.001 ~ "***",
      T ~ ""))
  
  est_all_1 <- tidy(model_1, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_2 <- tidy(model_2, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Spectators") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  est_all_3 <- tidy(model_3, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Luck",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  est_all_4 <- tidy(model_4, conf.int = TRUE, conf.level = 0.99) %>%
    mutate(model = "Reference Group: Routine Work",
           Group = "Workers") %>%
    filter(term %in% c("factor(treatment)1","factor(treatment)2"))
  
  
  all_treat_2 <- rbind(est_all_1, est_all_2, est_all_3, est_all_4) %>%
    dplyr::mutate(term = case_when(
      term == "factor(treatment)1" ~ "Routine Work",
      term == "factor(treatment)2" ~ "Complex Work"
    ))
  
  
  all_treat <- left_join(all_treat, all_treat_2, by = c("term", "estimate", "std.error", "statistic", "p.value", "model", "Group"))
  
  
  
  all_treat$term <- factor(
    all_treat$term, 
    levels = c("Complex Work","Routine Work"))
  
  
  
  all_plot <- all_treat %>%
    filter(model == "Reference Group: Luck") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+ 
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+ 
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") + 
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-40,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          axis.title.x=element_blank(),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.position = "none") +
    facet_wrap(.~model)
  
  
  
  
  
  reduced_plot <- all_treat %>%
    filter(model == "Reference Group: Routine Work") %>%
    ggplot(., aes(x=term, y=estimate, color = Group))  +
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.x, ymax=conf.high.x), width=0, size = 1.2)+
    geom_errorbar(position = position_dodge(width = 0.8), aes(ymin=conf.low.y, ymax=conf.high.y), width=0, size = 0.6)+
    geom_text(position = position_dodge(width = 0.8), aes(label=signif_stars), vjust = -0.35, size = 5, show.legend = FALSE)  +
    theme_bw()+ 
    geom_point(size=3, position = position_dodge(width = 0.8)) +
    labs(y = "Treatment Effect on Tax Rate",
         x = "") +
    scale_colour_grey() +
    geom_hline(yintercept=0, linetype = "dashed") +
    scale_y_continuous(breaks=seq(-40,10,10)) + 
    ylim(-40,10) +
    coord_flip() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.title=element_blank(),
          legend.background = element_blank(),
          legend.box.background = element_rect(colour = "black"),
          legend.position = "bottom")  +
    facet_wrap(.~model)
  
  
  combine_descriptive <- ggarrange(all_plot, reduced_plot, nrow = 2, heights=c(3,4))

  
####13. Perception Complexity Different Incomes####
  
  first_only <- vignettes_2 %>%
    mutate(first_decision = case_when(
      fl_7_do == "V1.$25K|V2.$50K|V3.$100K|V4.$500K" ~ b25,
      fl_7_do == "V1.$25K|V2.$50K|V4.$500K|V3.$100K" ~ b25,
      fl_7_do == "V1.$25K|V3.$100K|V2.$50K|V4.$500K" ~ b25,
      fl_7_do == "V1.$25K|V3.$100K|V4.$500K|V2.$50K" ~ b25,
      fl_7_do == "V1.$25K|V4.$500K|V3.$100K|V2.$50K" ~ b25,
      fl_7_do == "V1.$25K|V4.$500K|V2.$50K|V3.$100K" ~ b25,
    ))
  first_only$b25
  
  
  vignettes_test <- vignettes_2 %>%
    dplyr::select(b25,b50,b100,b500) %>%
    pivot_longer(cols = c(b25,b50,b100,b500)) %>%
    dplyr::mutate(name = case_when(
      name == "b25" ~ "25,000",
      name == "b50" ~ "50,000",
      name == "b100" ~ "100,000",
      name == "b500" ~ "500,000"
    ))
  
  
  predicted_rate <- vignettes_test %>%
    ggplot(., aes(x = factor(name), y=value)) +
    geom_point(position = position_dodge(width = 0.5),stat="summary",size = 3) + 
    geom_errorbar(position = position_dodge(width = 0.5),stat="summary", fun.data="mean_se", fun.args = list(mult = 1.96), width=0, size = 1.6) +
    geom_errorbar(position = position_dodge(width = 0.5),stat="summary", fun.data="mean_se", fun.args = list(mult = 2.56), width=0, size = 0.8) +
    theme_bw()+ 
    labs(y = "Predicted Complexity Perception",
         x = "Income Vignette") + 
    scale_color_grey() +
    theme(text = element_text(size=20),
          legend.text=element_text(size=15),
          plot.background = element_blank(),
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          legend.title=element_blank(),
          legend.background = element_blank(),
          legend.box.background = element_rect(colour = "black"),
          legend.position = c(.8,.9))
  

  vignettes_test %>%
    group_by(name) %>%
    summarise(mean_mean = mean(value))
  
  
  reorgan_data_vignettes <- vignettes_2 %>%
    dplyr::select(b25,b50,b100,b500,prolific_id) %>%
    pivot_longer(cols = -c(prolific_id))  
  summary(model_1 <- lm_robust( value ~ factor(name, levels = c("b25", "b50", "b100","b500")), reorgan_data_vignettes, clusters = prolific_id)) 
  
  
screenreg(l = list(model_1), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE)
  
####14. DAG####
  
  coord_dag <- list(
    x = c(c = 0, z = 1, y = 1, x = 2, p = 3),
    y = c(c = 0, z = -1, y = 1, x = 0, p = 0)
  )
  
  dagified <- dagify(y ~ c,
                     z ~ c,
                     x ~ y,
                     x ~ z,
                     p ~ x,
                     exposure = "c",
                     outcome = "p",
                     coords = coord_dag
  ) %>%
    tidy_dagitty(layout = "fr") %>%
    dag_label(c(y = "Effort Perception", z = "Skill Perception", 
                c = "Complex Work", x = "Deservingness/Fairness", p = "Tax Rate Preference")) 
  
  
  DAG <- ggdag_status(dagified, text = FALSE, 
                      use_labels = "label", edge_type = "link_arc", text_size =4, check_overlap = T) +
    theme_dag() +
    guides(fill = FALSE, color = FALSE)
  
  ####16. With covariates####
  
  spec_cov <- spec %>%
    mutate(college = case_when(
      d4 == "2-year college degree" ~1,
      d4 == "4-year college degree" ~1,
      d4 == "Doctoral degree" ~1,
      d4 == "Professional degree (JD, MD, MBA)" ~1,
      is.na(d4) ~ NA_real_,
      T ~0)) %>% 
    mutate(income = case_when(
      d5 == "Under $10,000" ~ 5000,
      d5 == "$10,000 - $20,000" ~ 15000,
      d5 == "$20,001 - $30,000" ~ 25000,
      d5 == "$30,001 - $40,000" ~ 35000,
      d5 == "$40,001 - $50,000" ~ 45000,
      d5 == "$50,001 - $60,000" ~ 55000,
      d5 == "$60,001 - $80,000" ~ 70000,
      d5 == "$80,001 - $100,000" ~ 90000,
      d5 == "$100,001 - $150,000" ~ 125000,
      d5 == "$150,001 - $200,000" ~ 175000,
      d5 == "$200,001 - $350,000" ~ 275000,
      d5 == "$350,001 - $500,000" ~ 425000,
      T ~ NA_real_)) %>%
    mutate(gender = case_when(
      d2 == "Male" ~ 1,
      d2 == "Female" ~ 0,
      T ~ NA_real_
    )) %>%
    mutate(party = case_when(
      d9 == "Republican Party" ~ 1,
      d9 == "Democratic Party" ~ 0,
      T ~ NA_real_
    )) %>%
    dplyr::select(d1, gender, income, d9, d6, d8_1, college, rd, treatment) %>%
    rename(age = d1,
           party = d9,
           employment = d6,
           left = d8_1) %>%
    mutate(treat = 1) %>%
    mutate(unemployed = case_when(
      employment == "Unemployed and looking for work" ~ 1,
      T ~ 0)) %>%
    mutate(retired = case_when(
      employment == "Retiree" ~ 1,
      T ~ 0)) %>%
    mutate(student = case_when(
      employment == "Student" ~ 1,
      T ~ 0)) 
  
  summary(model_1 <- lm(rd ~ factor(treatment) +gender + income + age + college + left + party + unemployed , spec_cov))
  
  summary(model_2 <- lm(rd ~ factor(treatment) +gender + income + age + college + left + party + unemployed , spec_red_cov))
  
  
  
  screenreg(l = list(model_1, model_2), include.ci = FALSE,
            include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
            booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 2, include.missings = FALSE, no.margin = TRUE, 
            label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
            custom.coef.map = list("factor(treatment)1" = "Routine Work",
                                   "factor(treatment)2"= "Complex Work",
                                   "gender" = "Male",
                                   "income"= "Income",
                                   "college"= "College Degree",
                                   "age"= "Age",
                                   "partyRepublican Party" = "Republican (Ref. Dem)",
                                   "partyOther" = "Other (Ref. Dem)",
                                   "partyDon't know" = "Don't Know (Ref. Dem)",
                                   "left"= "Left-Right",
                                   "unemployed" = "Unemployed"))
  
  
  
####15. Compare Workers and Spectators####
  #Merge datasets
  #Spectator, worker
  spec_cov <- spec_cov%>%
    filter(treatment == 0) %>%
    mutate(treat = 1)
  
  workers_cov <- workers%>%
    mutate(college = case_when(
      d4 == "2-year college degree" ~1,
      d4 == "4-year college degree" ~1,
      d4 == "Doctoral degree" ~1,
      d4 == "Professional degree (JD, MD, MBA)" ~1,
      is.na(d4) ~ NA_real_,
      T ~0)) %>%
    mutate(gender = case_when(
      d2 == "Male" ~ 1,
      d2 == "Female" ~ 0,
      T ~ NA_real_
    )) %>%
   mutate(income = case_when(
      d5 == "Under $10,000" ~ 5000,
      d5 == "$10,000 - $20,000" ~ 15000,
      d5 == "$20,001 - $30,000" ~ 25000,
      d5 == "$30,001 - $40,000" ~ 35000,
      d5 == "$40,001 - $50,000" ~ 45000,
      d5 == "$50,001 - $60,000" ~ 55000,
      d5 == "$60,001 - $80,000" ~ 70000,
      d5 == "$80,001 - $100,000" ~ 90000,
      d5 == "$100,001 - $150,000" ~ 125000,
      d5 == "$150,001 - $200,000" ~ 175000,
      d5 == "$200,001 - $350,000" ~ 275000,
      d5 == "$350,001 - $500,000" ~ 425000,
      T ~ NA_real_)) %>%
    dplyr::select(d1, gender, income, d9, d6, d8_1, college,rd, treatment) %>%
    rename(age = d1,
           party = d9,
           employment = d6,
           left = d8_1) %>%
    mutate(treat = 0) %>%
    mutate(unemployed = case_when(
      employment == "Unemployed and looking for work" ~ 1,
      T ~ 0)) %>%
    mutate(retired = case_when(
      employment == "Retiree" ~ 1,
      T ~ 0)) %>%
    mutate(student = case_when(
      employment == "Student" ~ 1,
      T ~ 0)) %>%
    mutate(treat = 0)
  
  
  both <- rbind(spec_cov,workers_cov)
  
  
  summary(model_balance <- lm(treat ~ gender + income + age + college + left + party + unemployed , data = both))
  
screenreg(l = list(model_balance), include.ci = FALSE,
         include.bic = TRUE, include.aic = TRUE, include.rsquared = T, include.variance = FALSE,  include.deviance = F,
         booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
         label = "table:main_treatment", caption = "Regression Results for Treatment Effects on Tax Rate for Top Earner", caption.above = TRUE, 
         custom.coef.map = list("gender" = "Male",
                                "income"= "Income",
                                "college"= "College Degree",
                                "age"= "Age",
                                "partyRepublican Party" = "Republican (Ref. Dem)",
                                "partyOther" = "Other (Ref. Dem)",
                                "partyDon't know" = "Don't Know (Ref. Dem)",
                                "left"= "Left-Right",
                                "unemployed" = "Unemployed"))
 